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Qh' Abstract: Based on the nonrelativistic QCD (NRQCD) factorization formalism, we 

fS ,■ investigate inclusive productions of two spin-triplet S-wave quarkonia pp — > 2J/ip + X, 

2T + X, and J/ip + T + X at the CERN Large Hadron Collider. The total production 

rates integrated over the rapidity (y) and transverse-momentum {pt) ranges \y\ < 2.4 and 

(S| ! p T < 50GeV are predicted to be a\pp -> 2J/^ + X] = 22 (35) nb, a\pp -)■ 2T + X] = 

24 (49) pb, and <r[pp — > J/ip + T + X] = 7 (13) pb at the center-of-momentum energy 

Q\ [ yfs = 7 (14) TeV. In order to provide predictions that can be useful in both small- and 

large-£>T regions, we do not employ the fragmentation approximation and we include the 

spin-triplet <S-wave color-singlet and color-octet channels for each quarkonium final state at 

^) • leading order in the strong coupling. The pt distributions of pp —> 2J/ip + X and 2T + X 

in the low-px region are dominated by the color-singlet contributions. At leading order in 

• • . the strong coupling, the color-singlet channel is absent for pp — > J/ip + T + X. Therefore, 

the process pp — > J/tp + T + X may provide a useful probe to the color-octet mechanism 

X ! of NRQCD. 
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1 Introduction 

Understanding the heavy-quarkonium production mechanism has been a longstanding 
problem in QCD, albeit considerable efforts both in theory and experiments. The color- 
singlet model [1-5] is most intuitive, and was the first attempt to describe the data in 
early days in QCD. This has been implemented as a leading-order (LO) contribution in the 
nonrelativistic QCD (NRQCD) factorization approach [? ] as vq — > 0, where vq is half 
the relative velocity of the heavy quark Q or heavy antiquark Q in the rest frame of the 
quarkonium H . Another phenomenological approach is the color-evaporation model [10? 
]. Although it is difficult to make the color-evaporation model systematic within QCD, 
this approach may capture an important point in the heavy-quarkonium production. This 



is especially the case, because of the general failure of local duality between quarks/gluons 
and hadrons. 

A salient feature of NRQCD is that color-octet contributions can be included in a sys- 
tematic fashion. The velocity-scaling rules (VSR) of NRQCD scale the numerical value for 
a long-distance NRQCD matrix element (0„ ( 2s+1 Lj)), which accounts for the transition 
probability of the heavy-quark-antiquark pair QQn{ 2s+1 Lj) to evolve into the quarkonium 
H in the asymptotic future, in powers of vq. Here, s, L, and J are the quantum numbers for 
the spin, the orbital angular momentum, and the total angular momentum of the QQ pair, 
respectively, and n = 1 (8) stands for the color-singlet (color-octet) state. In combination 
with the power counting in the strong coupling constant a s , physical observables can be 
expanded in double power series in a s and vq. In particular, in NRQCD the infrared (IR) 
divergence in the P-wave quarkonium \Qj decay in the color-singlet model is absorbed by 
renormalizing the color-octet matrix element (O g J ( 3 Si)), resulting in the decay rate free 
of an IR divergence [11]. 

The color-octet mechanism of NRQCD also enabled one to explain the large discrep- 
ancy between the data and the color-singlet-model prediction for prompt spin-triplet S- 
wave charmonium production at the Fermilab Tevatron [12] at large transverse momentum 
(pr)- It was possible by allowing a ccs( 3 5'i) pair, that was created from the fragmenta- 
tion of a virtual gluon, to make a transition into the spin-triplet S-wave quarkonium H. 
The corresponding matrix element (Og( 3 Si)} was determined from the large-py data. For 
the lower-pT region, other color-octet contributions ccs( 1 5'o) and ccs( 3 -Pj) as well as the 
color-singlet [cci( 3 Si)] channel may also contribute and a linear combination M r of the two 
matrix elements (Og( l So)} and (Og( 3 Pj)} was fit to the Tevatron data. The resolution 
of the large-surplus puzzle of prompt J/ip at large px by using the gluon fragmentation 
into ccg( 3 Si) was followed by the prediction that the polarization of prompt J/ip must be 
transversely polarized [13], which has not been confirmed by the data [14, 15]. 

However, the values for the color-octet matrix elements (Og ( 3 Si)), (Og ( 1 S'o)), and 
(Og ( 3 Pj)) are known not quite accurately [16] and there are several indications that these 
three matrix elements fit to the Tevatron data have been overestimated. Considerable 
higher-order corrections in a s to the color-singlet contribution to the prompt- J /ip produc- 
tion rate at the Tevatron have been reported although there is still a sizable room for the 
color-octet contribution, especially at large px [17, 18]. The nonobservation of strongly 
transverse polarization of prompt J/ip also indicates that the value for the matrix element 
(Og ( 3 Si)) has been overestimated. This puzzle remains unresolved in spite of recent the- 
oretical improvements [18, 20]. An analysis at next-to-leading-order (NLO) in a s of the 
inclusive-J/^ cross section at the B factories [21] showed that the color-singlet contribu- 
tion agrees with the data measured by the Belle Collaboration [23] within uncertainties. 
According to this analysis, the upper bound of the color-octet matrix elements (O g ( 1 5'o)) 
or (Og ( 3 -Po)) is much smaller than those determined by other experiments [24]. Recent 
analyses of inclusive J/ip photoproduction at NLO accuracies in a s show that the color- 
singlet contribution fails to describe various features of the data at HERA [25] and that 
the color-octet mechanism explains the HI data in spite of the poor knowledge of the 
color-octet matrix elements [27]. In summary, there are no conclusive constraints on the 



color-octet matrix elements, and it is still controversial if the color-octet mechanism makes 
a substantial contribution to the inclusive production of a quarkonium. It would be highly 
desirable to identify some processes which depend only on a very few color-octet matrix 
elements so that one can derive strong phenomenological constraints. 

Since the observation of exclusive double-quarkonium final states by the Belle Collabo- 
ration [28], production of double-quarkonium states in e + e~ collisions has lead remarkable 
progress in understanding the interplay of relativistic corrections and NLO corrections in 
a s [29? ? ? ? ]. The study has recently been extended to hadroproduction like the inclu- 
sive productions of double J/ip's, double T's 1 , and a Be* Be* pair at the Tevatron and the 
CERN Large Hadron Collider (LHC) [44-46]. For the double- J j tp production at LO in a s , 
the color-singlet contribution dominates over the color-octet one at pt ^ 8GeV [45]. On 
the contrary, in the case of the double-T production, the color-octet contribution dominates 
over the whole range of pr [45]. These predictions may be useful to study the color-octet 
mechanism in quarkonium production. In refs. [45, 46], the authors employed the gluon- 
fragmentation approximation in order to estimate the color-octet contribution. Although 
the approximation may give a reliable prediction at large pr, it should lose its predictive 
power at low values of px- 

In this work, we study inclusive productions of spin-triplet S-wave heavy-quarkonium 
pairs at the LHC in pp collisions at the center-of-momentum (CM) energies y/s = 7 and 
14TeV. The final states considered are double J/ip's, double T's, and J/tp + T. The 
analysis is carried out at LO in vq and in a s without employing the gluon-fragmentation 
approximation. This may allow us to provide a more reliable predictions for the production 
rates in the intermediate pt region, where the production rate is large in comparison with 
the large-pr region. The J/tp + T final state has a special feature that the color-singlet 
contribution is absent at LO in a s . Hence this process might be a clean probe to color-octet 
mechanism. Conversely, if J/tp + T events are not observed at the proposed level, it may 
lower the current upper bounds for the color-octet matrix elements significantly. 

This paper is organized as follows. In section 2, we describe basic strategies to compute 
the inclusive cross sections for pp — > 2 J/tp + X, 2T + X, and J/tp + T + X. We list various 
input parameters in section 3 and our predictions are given in section 4. We conclude in 
section 5 and provide relevant parton-level cross-section formulas in the appendix. 

2 Inclusive double-quarkonium production 

The NRQCD factorization formula for the differential cross section da of inclusive double- 
quarkonium production in proton-proton collisions, pp — > H\{pi) + #2(^2) + X, has the 
following schematic form: 

da[pp^H 1 + H 2 + X]= Y, fa/ P ®fb/ P ®da[ab^Q^ + Q n 2 *](0^)(0%), (2.1) 

a, 6, ni,n2 

where pi is the momentum of the quarkonium Hi, f a / p (x a , fj,) is the parton distribution 
function (PDF) for the parton a with the longitudinal momentum fraction x a with respect 



throughout this paper we suppress the identifier (IS) for T(IS'). 



Figure 1. Typical Fcynman diagrams for the nonfragmcntation contribution to pp — > 2H + X 
at LO in a s . Only gluon-gluon fusion diagrams are shown. The quark lines represent the charm 
(bottom) quark for H = J/ip (T). 




Figure 2. Typical Fcynman diagrams for the gluon- fragmentation contribution to doublc- 
quarkonium production processes pp — > 2 J/ip + X, 2T + X, and J/ip + T + X at LO in a s . 



to the proton, n is the factorization scale, the symbol <S> indicates the convolution over the 
partons' longitudinal momentum fractions x a and Xf,, and the summation is over all possible 
combinations of partons a and b. The parton-level differential cross section da[ab — > Q™ 1 + 
Q2 2 ] is the short-distance coefficient, which is perturbatively calculable in powers of a s . 
Here, Q™ 1 = QQ n { 2s+1 Lj) denotes the i-th QQ pair with the momentum pi and with the 
spectroscopic index m, which evolves asymptotically into H{. The dependence of the long- 
distance nature of the heavy quarkonium Hi is factored into the NRQCD matrix element 
(O^). The expression (2.1) is a power series in a s , v c , and Vf,. 

At LO in a s , only gg fusion and qq annihilation contribute to the parton processes 
ab — > Q™ 1 + Qg 2 ^ both Qi and Q2 are of the same flavor. The processes with the initial 
partons ab, ba = gq and gq are missing at this order. At the CM energies y/s = 7 and 
14 TeV of the LHC, one probes the small-x region of the PDF, where the gluon contribution 
dominates over the quark contents. Therefore, we ignore the parton processes initiated from 
qq, gq, and gq states and consider only the gg initial states in this work. 

2.1 Inclusive identical spin-triplet S'-wave quarkonium pair production 

Let us first consider inclusive identical spin-triplet 5-wave quarkonium pair production 
pp —7- 2H + X, where H is J/tp or T. In this case, the LO parton process is of order af . As 
is stated earlier in this section, we consider only the two-gluon initial states among various 
parton processes. The corresponding Feynman diagrams are classified into two groups: one 
is the nonfragmentation contribution, whose typical diagrams are shown in figure 1, and the 
other is the gluon- fragmentation contribution shown in figure 2. The sum of the two sets of 
diagrams can only make a gauge-invariant amplitude. In the limit that the invariant mass 
of Qi vanishes, the set of fragmentation diagrams approaches a gauge- invariant subset. 



According to VSR, the leading spectroscopic states of Q™ 1 and Q^ 2 in the velocity ex- 
pansion for H = J/ip or T production are both QQi( 3 S\) unless there is any enhancement 
factor for other spectroscopic states to compensate the power suppression compared to 
this color-singlet contribution. For the nonfragmentation diagrams in figure 1, the leading 
contribution must be the color-singlet channel QQi( 3 S\) + QQi( 3 Si). In the case of the 
gluon-fragmentation contribution QQs{ 3 S\) + QQg( 3 S±) in figure 2, the large kinematic en- 
hancement factor for the color-octet spin-triplet S- wave [QQs( 3 Si)] contribution, which is 
suppressed by a relative order v% compared to the color-singlet contribution, actually over- 
comes the suppression factor in the large-py region. Here, Q = c(b) for H = J/ip (T). How- 
ever, such enhancement factors do not appear in the mixed channels QQs( 3 Si) + QQi( 3 S±) 
and QQi( 3 Si)+QQs( 3 Si) that are suppressed by v% compared to the color-singlet channel. 
There are two other color-octet contributions, the spin-singlet S-wave [QQsC So)] and the 
spin-triplet P-wave [QQs( 3 Pj)], where J = 0, 1, and 2, which are suppressed by Vq and 
Vq, respectively, compared to the color-singlet contribution. Nevertheless, we ignore these 
two contributions because they do not have any large enhancement factor to compete with 
either QQi( 3 Si) or QQs( 3 Si) contribution. 

In summary, the contributions to the identical-quarkonium pair production that we 
consider in this work are QQi( 3 Si) + QQi( 3 Si) and QQ 8 ( 3 Si) + QQ 8 ( 3 Si). There are 31 
Feynman diagrams that contribute to the color-singlet channel gg — > QQi( 3 Si) + QQi( 3 Si) 
whose typical diagrams are shown in figure 1. For the color-octet channels gg — > QQs( 3 S\)+ 
QQs( 3 Si), there are 72 Feynman diagrams, some of which are shown in figures 1 and 2. 

The differential cross section in eq. (2.1) is applicable to both pp — > 2J/ip+X and pp — > 
2T + X if we substitute Q = c and b in the short-distance coefficients and H = J ftp and 
T in the NRQCD matrix elements, respectively. In the range where pt is not sufficiently 
large, the nonfragmentation contributions in figure 1 are not suppressed. Therefore, we 
compute the parton cross sections da[ab — > Q™ 1 + O^ 2 } for the relevant channels including 
all types of diagrams shown in figures 1 and 2. The parton-level cross sections for the 
color-singlet contributions are given in refs. [45, 46] and we have reproduced the results 
explicitly. The expression for the parton-level cross section for the color-octet channel 
gg — > QQs( 3 Si) + QQs( 3 Si) is given in appendix A. In the large-pr region, the color- 
octet gluon-fragmentation diagrams in figure 2 dominate and they can be computed by 
employing the fragmentation approximation of the gluon that fragments into a QQs( 3 Si) 
pair. In refs. [45, 46], such an approximation has been used to compute the cross section for 
the double-quarkonium production at large pt- In these references the parton cross sections 
for the real-gluon final states ab — > gg are convolved with the gluon-fragmentation function. 
In this work, we do not employ the fragmentation approximation and compute the complete 
set of order-a^ Feynman diagrams for QQi( 3 S\) + QQi( 3 Si) and QQs( 3 Si) + QQs{ 3 Si) 
channels. Therefore, our calculations provide predictions for these processes which can be 
compared with pt spectra in a wider range of pt including the small pt region. Further 
inclusion of the order-a^ color-octet channels that we have ignored in this work may improve 
the predictions in the intermediate-pT region, not modifying low- and large-py spectra 
significantly. This requires extensive calculations of formidably many Feynman diagrams 
and is beyond the scope of this work. 



(a) (b) (c) (d) 

Figure 3. Typical Feynman diagrams contributing to pp — > J/ip + T at LO in a s . 



2.2 J/ip + T production 

Next we consider the inclusive double-quarkonium production with different flavors, pp — > 
J /if) + T + X. Like the case of the same flavor, order-a 4 diagrams are leading in powers of 
a s . Typical Feynman diagrams for the gluon-initiated parton processes gg — > Q" 1 + Q^ 2 
are given in figure 3, where Q™ 1 (Q2 2 ) stands for the cc (bb) pair. At this order, the 
color-singlet channel cc~i( 3 Si) +bb~i( 3 Si) is absent. 

2.2.1 Color-octet contributions 

We first consider the color-octet contribution ccs( 3 Si) + bbsi^Si), whose "velocity-scaling 
factor" V 
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There are 36 Feynman diagrams that contribute to this channel and 
we show typical ones in figure 3. At large px the double- fragmentation contribution 
in figure 3 (a) dominates the cross section because of the large kinematic enhancement. 
Here, double fragmentation denotes that both cc and bb pairs are produced via gluon 
fragmentation. Figure 3(b) also represents mixed contributions ccs( 3 S±) + 6&i( 3 Si) and 
cci( 3 Si) + 66g ( 3 Si), each of which has 6 Feynman diagrams. Their velocity-scaling factors 
are v\ and v 4 , respectively, that are enhanced by either 1/u 4 or 1/u 4 compared to the chan- 
nel cc8( 3 5 , i) + 66s( 3 5i). Therefore, if pt is not large enough, then these mixed contributions 
must dominate over the color-octet channel cc8( 3 5i) + 66g( 3 S'i) by the enhancement factors 
l/i> 4 or 1/f 4 , while the double- fragmentation contribution dominates at large px- 

There are other contributions that are power suppressed compared to the contributions 
listed above. Such contributions are ccg( 3 Si) + 66§( 2s+1 Lj) and ccs( 2s+1 Lj) + 668( 3 Si), 
where 2s+l Lj = 1 Sq or 3 Pj. The corresponding Feynman diagrams representing these 
contributions are shown in figure 3(b), (c), and (d), with V = v^v^. Here, a, (3 = 3 for 
2s+1 Lj = 1 Sq and 4 for 3 Si and 3 Pj. The velocity-scaling factors of these contributions 
are approximately similar to that of the channel cc~g( 3 Si) + 66s( 3 5'i) but they do not have 
any double- fragmentation contribution. At small values of pt, these contributions are 
suppressed at least by either u 3 or v 3 compared to the mixed channels cci( 3 S\) + 668( 3 <Si) 
and ccs( 3 Si) + 661 ( 3 S\). Although the single-fragmentation channel in figures 3(b) and 
(c) may grow up at large px, that contribution is dominated by the double- fragmentation 
[figure 3 (a)] by a factor of (ruc/px) 4 or (jn^/pr) 4 . Hence, it is consistent to ignore ccs( 3 Si)+ 
bbs( 2s+1 Lj) and ccs( 2s+1 Lj) + 66s( 3 5i) channels over the whole pt range. 

The last color-octet contributions we can consider are the channels cc%( 2s+l Lj) + 
66g( 2s +1 L'j,), where 2s+1 Lj and 2s +l L'j, are 1 Sq or 3 Pj. A typical Feynman diagram of 
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Figure 4. Typical Fcynman diagrams for the color-singlct contribution to gg — > J/tp + T + X at 
order a®. 



these channels is shown in figure 3 (d). These contributions are, again, scaled by V = v^v^, 
where a,/3 = 3 for 2s+1 Lj = 1 So and 4 for 3 Pj. Thus these contributions are suppressed 
compared to the color-singlet one. Because they do not have any fragmentation contribu- 
tions, they are dominated by the double-fragmentation [figure 3 (a)] by (m c /pT) 4 (mft/pr) 4 
at large pr- Therefore, we ignore these channels, too. 

In summary, all of the channels that we consider for the process pp — > J/tp + T + X 
in this work are cc 8 ( 3 Si) + 66 8 ( 3 Si), cc 8 ( 3 5i) + ^i( 3 Si), and cci( 3 Si) + 66 8 ( 3 5i). The 
parton-level differential cross sections for these processes are given in appendix B. 

2.2.2 Color-singlet contributions 

As shown in figure 4, the LO color-singlet contribution to pp —> J/tp + T + X is of order 
a® and is suppressed compared to the color-octet contribution by a 2 s . The corresponding 
parton processes consist of 3 types: Type A, which is shown in figure 4(a), is gg — > 
cci( 3 Si) + bbi( 3 Si) + gg, where the initial gluons are attached to the charm- and bottom- 
quark lines one by one, a single virtual gluon connects the two quark lines, and each quark 
line emits a single real gluon. Type B is, again, gg — > cci( 3 Si) + bbi^Si) + gg, where 
both initial gluons are attached to a single quark line, the other quark line emits two real 
gluons, and a virtual gluon connects the two quark lines. A Feynman diagram for Type B 
is shown in figure 4 (b) and the two types A and B interfere at the amplitude level. Type 
C, which is shown in figure 4 (c), is gg — > cci( 3 Si) + 66i( 3 S*i), where the initial gluons are 
attached to the charm- and bottom-quark lines one by one and two virtual gluons connect 
the two quark lines simultaneously. 

If we consider the powers in a s only, then this color-singlet contribution must be 
suppressed compared to the color-octet processes listed in section 2.2.1. To assure that 
this argument is valid, we can make a rough estimate of the color-singlet contribution to 
the inclusive J/tp + T production rate at hadron colliders. 

We first classify the scaling behavior of the amplitude for the double-fragmentation 
diagram shown in figure 3(a). The propagator for the exchanged gluon has the scaling 
\jp\- The gluon propagators attached to the cc and bb pairs are of order l/mK,, and 
1/rriy, respectively. The scaling of the product of the two triple-gluon vertices must be 
the square of the typical momentum squared p\. Therefore, the resultant relative scaling 



of the double-fragmentation process is vt v b/( rn 'j/ib rn T:)^ wnere we have included the VSR 



suppression factor V = v\v\ 



,J b- 
In the case of Type A, there are four heavy-quark propagators and one gluon propagator 

which have large momentum transfer that account for the scaling factor l/p T . The phase 
space for the Qi + Q2 + 99 final state enhances the scaling factor by an order of p T in 
comparison with that for the two-body final state Q\ + Q2 of the double-fragmentation 
process. As a result, the scaling factor for Type A is l/p T - In a similar manner, we find 
that the scaling factor for Type B is the same as that of Type A. In Types A and B, there 
are two extra hard jets in the final states, and we have to include the coupling and the 
phase-space suppression factor ~ g 2 /(47r) 2 = a s /(4-/r) for each hard gluons. In addition, 
the color-octet contribution, which does not have hard jets, is distinguished from these 
color-singlet contributions that can simply be removed by imposing an appropriate veto. 
Type C diagrams involve finite box diagrams with two virtual gluons whose momentum 
must be of order the typical momentum transfer pr- A rough estimate of the scaling can 
be found by substituting a typical momentum transfer px to the two gluon propagators, 
four heavy-quark propagators, and the measure of the loop momentum. Furthermore, the 
one-loop amplitude has an additional suppression factor of g 2 /(47r) 2 = a s /(4ir). Therefore, 
the suppression factor for the color-singlet contribution relative to the color-octet double- 
fragmentation process is 

(4vr) 2 v£v* V PT Pt) ' { ' ' 

where we have included the strong-coupling suppression factor a 2 to the color-singlet con- 
tribution. Therefore one can argue that the color-singlet contributions to pp — > J/ip+Y+X 
will be suppressed enough compared to the color-octet ones. 

However, if pt is small, then the factors pt in this scaling should be of order ran and, 
therefore, the suppression factors for the color-singlet contribution relative to the mixed 
contributions cci( 3 Si) + bb^S\) and ccs( 3 Si) + 6&i( 3 Si) are approximately estimated to 
be ~ a 2 /[(47r) 2 ^] or a 2 /[(47r) 2 i^] in the small px region. These factors are much less than 
order 1 if we assume t> 2 ~ 0.3 and u 2 ~ 0.1 and we use the renormalization scale of order 
mft . Thus we may conclude that the color-singlet contribution is suppressed compared to 
the color-octet contribution coming from the mixed channels even in the small pt region. 
This rough estimate might fail as pt approaches to zero. Then, our prediction for the 
ccs( 3 'Si)+bbs( 3 'S\) contribution may be contaminated by order-a^ color-singlet contribution 
at small px- The problem can be resolved by introducing a lower px cut pt > 5 GeV, 
where the double fragmentation rises up. The suppression of the color-singlet contribution 
to pp —?• J /ip + T + X even in moderate values of px has not yet been predicted, because 
previous analyses with the fragmentation approximation are valid only for the large-py 
region. In addition, the color-singlet contribution is easily distinguishable by imposing 
a veto that the final state must not include hard jets. This is the main motivation for 
investigating pp — > J/ijj + T + X in this work as a clean probe to the color-octet mechanism 
at the LHC. This point has been recently reported in ref. [47]. 



2.3 Gluon-fragmentation approximation 

In inclusive single-hadron production in pp collisions, a further factorization happens if 
Pt is sufficiently large [48]. In the case of the 3 Si heavy-quarkonium production, the 
fragmentation of a gluon into a QQsi^Si) pair dominates the production rate at large 
Pt [12]. If the two quarkonia produced in pp —> Hi + Hi + X both have large transverse 
momenta, then one can guess that such a factorization can be generalized so that one might 
be able to write the inclusive double-quarkonium production cross section in the form: 

d<Ji \pp -)■ H x + H 2 + X] = f g/p (g> f g/p <g) da[gg -> gg] <g) D g ^ Hl <g) D g ^ H2 , (2.3) 

where the subscript in dai indicates the fragmentation approximation, da[gg — > gg] is the 
parton-level cross section for gg — > gg, and D g ^H(z, flf) is the fragmentation function 
for a gluon to fragment into a QQsi^Si) pair that evolves asymptotically into the heavy 
quarkonium H = J/tj) or T. Here, z is the longitudinal momentum fraction of H relative to 
the fragmenting gluon and /if is the factorization scale for the fragmentation. An educated 
guess is that the factorization scale /if j for the fragmentation of a gluon into the quarkonium 
Hi can be chosen to be of order the transverse mass niTi = (wif + PxiY ■> where rrii and 
PTi are the mass and the transverse momentum of Hi, respectively. 

At the threshold /if = 2niQ, D g ^H{z,m = 2?tiq) is calculable perturbatively after 
factoring out the NRQCD matrix element (Og ( 3 <Si)}. The fragmentation function of a 
gluon that fragments into the 3 Si quarkonium via QQs( 3 Si) pair is known up to NLO in 
a s at the threshold [49, 50]. The expression at LO in a s is given by [? ] 



7ra g 
24m^, 



D g ^ H (z,2m Q ) = — U(l ~ z)(Of ( 3 5f)). (2.4) 



In order to take into account multiple emissions of collinear gluons before creating a QQ 
pair, it is necessary to evaluate the fragmentation function at the factorization scale /if ~ 
Pt 3> rriQ. This step can be carried out by making use of the Altarelli-Parisi evolution 
equation [51-53]. 

Very recently, Qiao, Sun, and Sun [46] considered the inclusive process pp — > 2J/ip + X 
and Li, Zhang, and Chao [45] carried out extensive studies of various double-quarkonium 
production processes that include pp — > 2J/ip + X and 2T + X, where the gluon-gluon 
initial states were mainly considered as are in this work. While the color-single contribution 
was computed by using the same way that we employ here, the authors of these papers 
used the gluon-fragmentation approximation to estimate the color-octet contributions to 
pp —?• 2J/tp + X [45, 46] and 2T + X [45] 2 . However, in refs. [45, 46] the fragmentation 
function D 9 ^h(z, /Xf) was not evaluated at /tf ~ pt but at /if = 2itiq. In fact, in order to 
make the fragmentation approximation valid, the fragmenting gluon has energies of order 
Pt 3> rriQ where the strong coupling is significantly smaller. Therefore, one may guess 



2 Note that our analysis on the process pp — > J/ip + T + X is the first NRQCD-based study including 
both color-singlet and -octet channels that does not employ the gluon-fragmentation approximation, while 
the authors of ref. [10] also have commented about the process at low-energy hadron collisions by making 
use of the color-evaporation model. 



that the predictions given in refs. [45, 46] must have been overestimated. On the other 
hand, once we employ the gluon-fragmentation approximation to compute the color-octet 
contributions, the prediction cannot be extended to a lower pt region, where the octet 
contribution acquires IR divergences. At small pt, the fragmentation contribution loses its 
predictive power. Therefore, the predictions in refs. [45, 46] do not extend to lower values 
of pt- Thus, our choice of not employing the gluon-fragmentation approximation in this 
work enables us to predict the pt spectrum to a wider range whose lower bound extends 
to pt = 0. 

2.4 QED contribution 

As another color-singlet contribution to the double-quarkonium production at hadron col- 
liders, one may consider the parton-level process qq —> QQi( 3 Si) + QQi( 3 S±) . This subpro- 
cess may acquire a large kinematic enhancement factor due to the double photon fragmen- 
tation, qq — > 7*7* followed by 7* — > QQi^Si) for each virtual photon. As the double- J/ip 
production in e + e~ annihilation into two virtual photons, this process has additional loga- 
rithmic enhancement in the forward region [? ]. However, in hadroproduction, it is very dif- 
ficult to analyze that kinematic region. In addition, the rough estimate of this QED contri- 
bution to the double-gluon-fragmentation contribution is ~ e^e^a 4 / (a 4 t> 4 i; 4 ) ~ O(10~ 5 ) 
with the electromagnetic coupling a and electric charges e„, e c , and e& of the initial, charm, 
and bottom quark, respectively, which is negligible numerically. This (/(/-initiated QED con- 
tribution is further suppressed compared to the gluon-initiated processes. Therefore, we 
do not consider the QED contribution in this work. 

3 Numerical Analysis 

Based on the formalism described in section 2, we are ready to carry out the numerical 
calculation of the production rates for the double-quarkonium production processes pp — > 
2J/tp + X, 2T + X, and J/ip + T + X at the LHC with the CM energies y/s = 7 and 14TeV. 
In this section, we describe our choice of various input parameters and the strategies of the 
numerical evaluation of these cross sections. 

3.1 Parameters involving NRQCD factorization 

In order to evaluate the production rates in eq. (2.1), we have to know the values for the 
long-distance NRQCD matrix elements (0^( 3 Si)} for H = J/ip and T, where n = 1 or 8. 
The color-singlet matrix element (Oi( 3 Si)) h for the spin-triplet S-wave heavy quarkonium 
decay is usually determined from the leptonic decay rate of H, which is the most precisely 
measured value involving (O n ( 3 Si))# 3 . Under the vacuum-saturation approximation, the 
decay matrix element is approximately a third of the production matrix element up to 



3 There has been a great progress in precise determination of the color-singlet matrix element of the 
J/ip meson due to the introduction of a new technique to deal with the relativistic corrections and its 
resummation in conjunction with the one-loop QCD correction to the electromagnetic decay rate [54-58]. 
An analogous method has been used to determine the matrix elements for corresponding bottomonium 
states [59]. 
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corrections of order v% [? ]: (0^( 3 Si)} = 3(Oi( 3 »Si)}# + 0(vq), where the factor of 3 
on the right side of the equality stands for the spin-multiplicity factor (2 J + 1) for the S- 
wave spin-triplet state. We quote the following values for the color-singlet NRQCD matrix 
elements for J/ip and T: 

(0^( 3 5i)> = 1.3 GeV 3 [56], (3.1a) 

(0?"( 3 5i)) = 9.2 GeV 3 [59]. (3.1b) 

The color-octet matrix element (O g ( 3 Si)) has been fit to the px spectrum of prompt J/ip 
production rate at the Tevatron in the large-p^ region [13]. The matrix element (Oj( 3 Si)) 
has also been fit [61] to the Tevatron data and used for the polarization analysis [62]. 
Various determinations of these matrix elements can be found, for example, in ref. [63]. 
We employ the following values for the 3 <Si color-octet NRQCD matrix elements for J/ip 
and T: 

(Og /V, ( 3 5i)) = 3.9 x 10" 3 GeV 3 [13], (3.2a) 

(Oj( 3 Si)) = 1.5 x KT 1 GeV 3 [63]. (3.2b) 

The short-distance coefficient, which is the parton-level differential cross section da in 
eq. (2.1), depends on the heavy-quark mass mq and the strong coupling a s . For mq, we 
take ra c = 1.5 GeV for the charm quark and m& = 4.7 GeV for the bottom quark. In the 
NRQCD factorization formula, the short-distance coefficients are expressed in terms of raq 
instead of the meson mass mjj- This may give rise to relativistic corrections due to the 
difference between rau and 2raq at higher orders in vq. However, because we carry out 
our calculation at LO in vq, we can ignore such corrections and put ran = Zmq in our 
analysis. 

3.2 Other input parameters 

In addition to the NRQCD factorization, the cross section formula (2.1) involves the fac- 
torization of the long-distance PDF's, f a / p and fb/ p , and the short-distance parton-level 
cross section da with the factorization scale fj,. For the scale //, we take the transverse 
mass n = rax = (4m?, + p T ) 1 ' 2 , which is a conventional choice 4 . In general, the transverse 
momenta for the two heavy quarkonia in the final state can be different. However, at LO 
in a s , there are no additional hard jets and the two final-state quarkonium pairs have the 
same px and, therefore, px and rax are defined unambiguously. As a specific choice of the 
PDF in eq. (2.1), we employ the CTEQ6L parametrization [64]. 

In evaluating a s , we set the renormalization scale to be rax so that a s = a s (/j, = rax)- 
In order to make our numerical evaluation of a s consistent with the CTEQ6L parametriza- 
tion, we use the NLO formula for the running coupling constant a s (/x) by setting a s (fi = 
M z ) = 0.118 with A 4 = 326 MeV [64]. 

The kinematic region that we study in this work covers the rapidity range |y| < 2.4. 
Our results for the color-octet contributions are finite as well as the color-singlet contribu- 
tion at px = 0, since we do not adopt the gluon-fragmentation approximation. Therefore, 



4 See, for example, ref. [13]. 
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Figure 5. The differential cross sections for pp — ► 2J/V> + 1 at (a) \/s = 7TeV and (b) 14TeV 
in units of nb/GeV as functions of pt integrated over the rapidity range \y\ < 2.4. The solid, 
dashed, and dotted curves represent the total, color-singlet [cci( 3 5*i) + cci( 3 Si)], and color-octet 
[ccs( 3 Si) + ccg( 3 Si)] contributions, respectively. 



we set the lower bound of px to be that is distinguished from the previous studies in 
refs. [45, 46], where they considered the range 3GeV < pr < 50GeV. 

4 Predictions 

In this section, we provide our predictions of the differential cross sections for the double- 
quarkonium productions pp —> 2J/tp + X, 2T + X, and J/tp + T + X at CM energies y/s = 
7 and 14 TeV. Detailed strategies and the values for the input parameters for the numerical 
analysis have been given in sections 2 and 3, respectively. 

4.1 da/dp T for pp -> 2J/i/j + X and 2T + X 

Our predictions for the pt dependence of the differential cross sections da/dpxlpp — >• 2J/ip+ 
X] at y/s = 7 and 14 TeV integrated over the rapidity range \y\ < 2.4 are shown as solid 
curves in figures 5 (a) and (b), respectively, and those for da/dprlpp ~^ 2~F + X] are shown 
in figure 6. As is stated in section 2.1, we have considered only order-o^ gluon-initiated 
parton processes: the color-singlet contribution through gg — > QQi(^Si) + QQ\( 3 Si) and 
the color-octet contribution gg — > QQsi^Si) + QQs( 3 Si) for these identical-quarkonium 
pair productions. 

As far as the shape is concerned, the px spectra at the two CM energies \/s = 7 
and 14 TeV are essentially the same. The ratio of the value at yfs = 14 TeV to that 
at y^ = 7 TeV is about a factor of 1.5 (2) for pp — > 23 '/tp + X (2T + X) near the peak, 
where the color-singlet contribution dominates. As pt increases, the ratio increases and the 
gluon- fragmentation contribution dominates. The pt spectra in figures 5 and 6 show that 
the color-singlet channel dominates at small px while the color-octet channel dominates at 
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Figure 6. The differential cross sections for pp — > 2T + X at (a) ^/s = 7 TeV and (b) 14 TeV 
in units of pb/GcV as functions of pr integrated over the rapidity range |y| < 2.4. The solid, 
dashed, and dotted curves represent the total, color-singlet [bbi( 3 Si) + 66i( 3 5i)], and color-octet 
[bbs( 3 Si) + bbs^Si)} contributions, respectively. 



large pt- The crossovers are placed at px ~ 16 GeV for pp — > 2J/ip-\-X and at px ~ 24 GeV 
for pp — > 2T + X at y/s = 7 and 14 TeV both. Another noticeable feature of the spectra is 
that the values at pt = vanish: 



lim - — [pp 
pt-H) dp? 



2H + X] = 0. 



(4.1) 



The reason is that, at order a^v%, all of the diagrams we consider in this work do not 
contain any real gluons in the final state and both color-singlet and -octet contributions 
are free of IR divergences 5 . In the low-px region, the differential cross sections increase 
as pt increases until they reach the maximum values da/dpxipp —> 2J/ip + X] = 11.2 
(17.3) nb/GeV at p T = 1.1 GeV and da/dp T [pp ->■ 2T + X] = 4.1 (8.0) pb/GeV at p T = 3 
GeV for ^ = 7 (14) TeV. 



y/s\a (nb) 


CC^S^ + CC^S!) 


ccseS^+ccs^S^ 


total 


7 TeV 

14 TeV 


22.3 

34.8 


0.011 
0.019 


22.3 

34.8 



Table 1. The total cross sections <r[pp —> 2J/ip + X] at yfs = 7 and 14 TeV integrated over the 
ranges \y\ < 2.4 and pt < 50 GeV in units of nb. The three columns represent the color-singlet, 
color-octet, and total contributions. 



5 However, if there is at least one real gluon in the final state, then in the limit that the final gluons have 
vanishing momenta, IR divergences may appear in the color-octet contributions even at LO Vq. In the case 
of the color-singlet contributions such an IR divergence may grow at higher orders in vq even at LO in a a ■ 
For more discussions, we refer the readers to refs. [65-68]. 
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Vs\<? (pb) 


bb l ( 3 S 1 ) + bb 1 ( 3 S l ) 


66 8 ( 3 5 1 ) + 66 8 ( 3 5 1 ) 


total 


7TeV 
14TeV 


24.1 
47.9 


0.27 
0.60 


24.4 
48.5 



Table 2. The total cross sections u[pp — > 2T + X] at y/s = 7 and 14TeV integrated over the 
ranges \y\ < 2.4 and pr < 50GeV in units of pb. The three columns represent the color-singlet, 
color-octet, and total contributions. 

Next we compute the total cross sections a\pp — > 2J/ip + X] and a\pp — > 2T + X] by 
integrating eq. (2.1) over the ranges \y\ < 2.4 and pt < 50GeV. As shown in tables 1 and 
2, they are a[pp ->■ 2 J/V> + X] = 22 (35) nb and cr[pp -> 2T + X] = 24 (49) pb at y/s = 7 
(14) TeV. Because the differential cross section is dominated by the color-singlet channel 
near the peak, the total cross section is dominated by the color-singlet contribution. 
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Figure 7. The differential cross sections for pp —$■ J/ip + T + X at (a) s/s = 7 TeV and (b) 14 TeV in 
units of pb/GeV as functions of px integrated over the rapidity range \y\ < 2.4. The solid, dashed, 
dashed-dotted, and dotted curves represent the total, cci( 3 Si) +bbs( 3 Si), ccs( 3 Si) + bbi( 3 Si), and 
ccs( 3 Si) + bbs( 3 Si) contributions, respectively. 



4.2 da/dpT for pp -> J/ip + T + X 

The predictions for da /dp? [pp — > J/^p + T + X] at ^fs = 7 TeV and 14 TeV integrated over 
the rapidity range \y\ < 2.4 are shown as solid curves in figures 7 (a) and (b), respectively. 
In this process, we consider the following order-o^ contributions: the color-octet channel 
gg — > ccg( 3 5i) + 66g( 3 S'i) and the two mixed channels gg — > cci( 3 S\) + bbs( 3 Si) and 
gg — > cc8( 3 £i) + bb\( 3 Si). The reason for neglecting the color-singlet contribution and the 
remaining color-octet channels is given in section 2.2. 

The pt spectrum for pp — > J/ip+T+X is similar to those of pp — > 2 J/ip+X and 2T+X. 
The differential cross section vanishes at pt = and rapidly increases until it reaches the 
maximum value da/dpT[pp — > J/ip + T + X] = 1.6 (2.9) pb/GeV at pr = 1.6 (1.5) GeV 
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for y/s = 7 (14) TeV. Then it monotonically decays as pt increases. Near px = 0, the 
cci( 3 Si) + 66g( 3 S'i) channel dominates and that is the largest contribution for px < 4 GeV. 
At large px, the color-octet channel ccs( 3 5i) + bbs^Si) dominates and the channel is the 
largest contribution for px > 6 GeV. In the remaining region, about 4 GeV < pt < 6 GeV, 
the color-octet channel ccs( 3 S'i) + bb^S\) and the two mixed channels cci( 3 Si) + bb$( 3 Si) 
and cc%( 3 Si) + &6i( 3 Si) compete together. 



^W P b) 



cci( 3 5i) + 66 8 ( 3 5i) cc 8 ( 3 5i) + 66i( 3 5i) cc&( 3 Si) + && 8 ( 3 -Si) total 



7 TeV 

14 TeV 



3.18 1.95 1.63 6.76 

6.00 3.72 3.36 13.08 



Table 3. The total cross sections a[pp — > J/tp + T + X] at \/i = 7 and 14 TeV integrated over the 
ranges \y \ < 2.4andpr < 50 GeV in units of pb. The four columns represent the cci( 3 Si) + bb$ ( 3 <Si), 
cc 8 ( 3 Si) + 66i( 3 5i), cc 8 ( 3 <Si) + bb$( 3 Si), and total contributions. 

We compute the total cross section a\pp — )• J/tp + T + X] by integrating eq. (2.1) 
over the ranges \y\ < 2.4 and pt < 50 GeV. As shown in table 3, the total cross section 
at y^i = 7(14) TeV is a[pp -» J/tp + T + X] = 6.8 (13) pb. The total cross section for 
pp — > J/ip + T + X without any px cut is smaller than those for 2J/ip and 2T final 
states by factors of ~ 3000 and ~ 4, respectively. We notice that most of the rates for 
the 2J/i/j and 2T cases are concentrated in the low px regions, where the color-singlet 
contribution dominates. In order to probe the color-octet contribution more accurately, 
we had better impose a lower pt cut to remove most of the color-singlet contribution. The 
cut for the 2 J/tp or 2T case should be around the place where the QQs( 3 Si) + QQs{ 3 Si) 
contribution overtakes the QQi( 3 Si) + QQi( 3 S\) contribution. According to figures 5 and 
6, the cuts pt > 16 and 24 GeV are appropriate for 2J/tp and 2T final states, respectively. 
In those cases, a\pp —?■ 2 J/tp + V]| PT > 16GcV = 0.09 (0.2) pb, a[pp -> 2T + V]| pT > 24 GeV = 
0.02(0.05) pb, respectively, at %fs = 7 (14) TeV. As is described in section 2.2, the whole 
cross section a[pp — > J /tp + T+X] = 6.8 (13) pb depends on the color-octet matrix elements 
and there is no need to impose an additional pt cut to remove the color-singlet contribution. 
Therefore, we conclude that the pp — > J/tp + T + X channel is the most sensitive to the 
color-octet matrix elements among the three double-quarkonium final states. 

4.3 Comparison with previous calculations 

We compare our predictions in sections 4.1 and 4.2 with previous calculations [45, 46] 
based on the gluon-fragmentation approximation. We find that our prediction is distin- 
guished from the previous calculations in the shape of the pt spectrum and these previous 
calculations severely overestimate the rates all over the range of px- 

First of all, while our pt spectra in figures 5-7 have maximum values around pt ~ 2, 
4-5, and 2-3 GeV for pp — > 2J/tp+X, 2T + V, and J/tp + T + X, respectively, the previous 
predictions in refs. [45, 46] do not have these peaks. The values of pt at which the peaks 
appear are of order mjj^ for pp — > 2 J/tp + X and pp — > J/tp + T + X and of order m-j for 
pp —7- 2T + X. This is the region where the fusion contributions are important. Neglecting 
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fusion diagrams in the fragmentation approximation cannot reproduce the correct shape 
of pt spectrum at transverse momenta less than a few run. 

In addition, those previous predictions based on the gluon-fragmentation approxima- 
tion overestimate the rates all over the pt range. At small pr, the fragmentation approxi- 
mation must break down and, therefore, the approximation severely overestimate the rate 
in this region. The approximation actually diverges as px — > while our prediction vanishes 
in that limit. In fact, the fragmentation approximation must be a good approximation if 
Pt is sufficiently large. At large pr, the predictions in refs. [45, 46], again, severely overes- 
timate the rate. The reason is that in these references the authors fixed the factorization 
scale for the gluon fragmentation to be fj,{ = 2rriQ, which leads to large cross sections in 
comparison with our predictions given in figures 5 and 6 in which we use conventional 
choice /if = tut- 

We can compare our results for the double-quarkonium production of the same fla- 
vor, pp — > 2J/ip + X and pp — > 2T + X with those for the single-quarkonium produc- 
tion associated with a heavy-quark pair with the same flavor, pp — > J/ip + cc + X and 
pp — > T + bb + X [69]. In both cases, the color-octet contributions are of the same order in 
a s as those for the color-singlet channels. Therefore, if there are no kinematic enhancement 
factors such as fragmentation, then the color-octet contributions must be suppressed. The 
results shown in figures 5 and 6 are consistent with this expectation that the color-singlet 
contribution dominates over the color-octet contribution unless pt is extremely large. Com- 
paring our results for the production rate with those in ref. [69] which contains only the 
color-singlet contribution, we find that the production rate for the single-quarkonium pro- 
duction associated with a heavy-quark pair with the same flavor at the LHC is by about 
a factor of 10 3 larger than that for the double-quarkonium production at the LHC, where 
we have taken into account the decay of each heavy-quarkonium into a muon pair. Such 
a suppression should be due to the restrictions in the phase space and the bound-state 
formation. However, the double-quarkonium production has clear signals, while the single- 
quarkonium production associated with an open heavy-quark pair may suffer from large 
backgrounds. 

5 Discussions 

Based on the NRQCD factorization formalism, we have obtained the px spectra and the 
total cross sections for pp — > 2J/ip + X, 2T + X, and J/ip + T + X at CM energies yfs = 
7 TeV and 14 TeV. The total production rates integrated over the rapidity and transverse- 
momentum ranges \y\ < 2.4 and pt < 50 GeV are predicted to be a\pp — > 2J/ip + X] = 22 
(35) nb, a[pp ->■ 2T + X] = 24 (49) pb, and a[pp ->■ J/i/> + T + X] = 7 (13) pb at ^/s = 7 
(14) TeV. 

For the double-quarkonium production of the same flavor, pp — > 2J '/tp + X and 
2T + X, production rates are dominated by the color-singlet contribution at low px- 
In order to probe the color-octet contribution, which dominates at large pt, more ac- 
curately we have imposed lower px cuts where the color-octet contribution overtakes 
the color-singlet one. The results are a[pp — > 2J/ip + AT]|p T >i6GeV = 0.09 (0.2) pb and 
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a[pp — > 2T + X]| pT >24GeV = 0-02 (0.05) pb. Assuming the integrated luminosity ~ lOOfb -1 
at y/s =14TeV and considering the branching fractions B[J/ip —> /i + /x~] = 5.93% and 
B[T — > /^ + /U~] = 2.48% [70], we expect that approximately 70 (3) double- J ftp (-T) events 
can be observed by tagging muon pairs under the cuts pt > 16 (24) GeV for pp — > 2J/ip+X 
(2T + X). The results indicate that it is very difficult to study the color-octet mechanism 
by making use of the processes pp — > 2 J '/ip + X and pp — > 2T + X at the LHC. 

In the case of pp — > J/ip + T + X, the color-singlet channel is absent at LO in a s , 
and is highly suppressed relative to the color-octet contributions [ see eq. (2) ]. The two 
channels cci( 3 Si) + bt>s( 3 Si) and cc$( 3 Si) + bbs( 3 Si) dominate the production rate at small 
and large values of px, respectively. The contribution ccs( 3 S'i) + 6&i( 3 Si) is comparable 
to those two listed above only around the region 4 GeV < pr < 6 GeV. Assuming the 
integrated luminosity ~ 100 fb~ at yfs =14TeV and considering the branching fractions 
B[J/ip -> /i + //~] = 5.93% and B[T -> /U+// _ ] = 2.48% [70], we expect that approximately 
1900, 520, and 160 events can be observed by tagging muon pairs under the cuts px > 
0, 5, and 10 GeV, respectively, at the LHC. Improving the acceptances for J/ip and T by 
extensive Monte Carlo studies of final-muon pairs, one may expect that the observation 
of the events and the determination of the matrix elements can be quite promising in the 
near future. The number of events may be increased by a factor of 4 if one includes the 
e + e~ decay modes of J/ip and T. It may also be improved by including the subprocess 
qq — > Hi + H2 via two-gluon exchanges that is neglected in this work. 

It is well known that, in inclusive single-quarkonium production in hadron collisions, 
there are large corrections at NLO in a s . For example, the NLO corrections to the color- 
singlet contribution to the inclusive J/ip production enhance the rate by an order of mag- 
nitude [17, 18, 71] especially at large pr- Therefore, one can worry that NLO corrections 
in a s may spoil the LO prediction for the pt spectra presented in this work. In inclusive 
single-quarkonium production, NLO subprocesses include gg — > J/ip + gg, gg — > J/ip + qq, 
9Q(q) ~^ J/tP + Q9(q), an d gg — > J/tp + cc. Although there is a large enhancement from 
i-channel gluon-exchange diagrams, it is also true that each of new diagrams has signifi- 
cant contributions, piling up the corrections to modify the rate by an order of magnitude. 
In the case of double-quarkonium production of the same flavor, new NLO subprocesses 
include t-channel gluon-exchange diagrams such as gg — > 2H + g and gq(q) — > 2H + q(q), 
where H = J/ip or T. But the number of new channels are quite limited compared to 
the single-quarkonium production case. Therefore, we expect that the NLO corrections to 
the color-singlet contribution to the double-quarkonium production of the same flavor may 
enhance the rate significantly, but not as dramatically as the single-quarkonium produc- 
tion, while observation of the color-octet mechanism might be significantly affected. This 
argument should be tested by a forthcoming quantitative analysis of the NLO corrections 
to the color-singlet and color-octet contributions. 

In the case of J/ip + T production, the suppression factor of the color-singlet channel 
to the color-octet one is a^/p^p which is significantly smaller than the corresponding fac- 
tor l/p\- for the single-quarkonium production, in which the color-singlet and color-octet 
contributions are of the same order (a 3 ) at LO. We anticipate that our arguments can be 
tested quantitatively by measurements and by explicit calculations of complete order-a^ 
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contributions to the color-singlet channel in near future. 
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A Parton-level cross section for gg — ► QQs( 3 Sx) + QQs( 3 Sx) 

In this appendix, we present the parton-level differential cross sections for pp — > 2J/ip + X 
and pp — > 2T + X. We first define the Mandelstam variables for the parton-level process 
g(kx)g(k 2 ) — > QQ ni (px) + Q'Q' n J(P2)-, where the momentum for each state is given in the 
parentheses following the state and QQ ni and Q'Q' n2 evolve into the quarkonium states H\ 
and H 2 , respectively. Each state is on its mass shell so that k\ = k\ = 0, p\ = (2mg) 2 , 
and p\ = (2toq') 2 . At LO in vq and vq>, the meson masses are set to be m# x = 2niQ and 

2itiq. The Mandelstam variables s, i, and 



{pi+p 2 f, (A.la) 

(k 2 -p 2 ) 2 , (A.lb) 

(k 2 -pi) 2 - (A.lc) 

In pp —)■ 2H + X for H = J/ijj or T we consider the two channels gg — > QQii^Si) + 
QQi^Sr) and gg -> QQ 8 ( 3 ^i) + QQ 8 ( 3 Si), where Q = c (b) for H = Jfij, (T). The 
parton-level cross sections for gg — > QQi(^Si) + QQi(^Si) can be found in refs. [45, 46]. 
That for the subprocess gg — > QQs^Si) + QQs( 3 Si) is given by 

f to - gACft) + gAC*)] = ^-^^— w p>^, (a. 2) 



ITIH2 = 2m,Qi . If H\ - 


= H 2 , then niH 1 = rnn 2 


u are defined by 






s = (h + k 2 f 




t = (kx - px? 




u = (kx -p 2 ) 2 
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s 



where mu = 2toq and the coefficients Oj's are given by 

a = 243s 4 * 2 (s + i) 2 (s 2 + si + ff, 

01 = -162s 3 P(s + *)(s 2 + at + t 2 )(9s 5 + 19s 4 i + 20s 3 * 2 + 7s 2 * 3 - 3si 4 - 6i 5 ), 

a 2 = t(243s n + 3951s 10 * + 6714s 9 * 2 + 14420s 8 * 3 + 179582s 7 * 4 + 919446s 6 * 5 

+2488136s 5 * 6 + 4132862s 4 * 7 + 4395900s 3 * 8 + 2933988s 2 * 9 + 1119744s* 10 

+186624* 11 ), 
a 3 = -2*(57s 10 + 1233s 9 * + 46541s 8 * 2 + 513120s 7 * 3 + 2646793s 6 * 4 + 7942109s 5 * 5 

+15041136S 4 * 6 + 18324922s 3 * 7 + 13942380s 2 * 8 + 6013080s* 9 + 1119744* 10 
a 4 = 2(935s 10 + 9398s 9 * + 117747s 8 * 2 + 1103652s 7 * 3 + 6182220s 6 * 4 + 21423546s 5 * 

+47491450S 4 * 6 + 67574132s 3 * 7 + 59508939s 2 * 8 + 29339460s* 9 + 6158916* 10 ), 
a 5 = -2(8039s 9 + 112887s 8 * + 1157014s 7 * 2 + 7632256s 6 * 3 + 31876569s 5 * 4 

+85147430s 4 * 5 + 144700858s 3 * 6 + 150182520s 2 * 7 + 85844772s* 8 + 20531880* 9 ), 
a 6 = 2(43072s 8 + 638490s 7 * + 5393635s 6 * 2 + 28486982s 5 * 3 + 94986651s 4 * 4 

+198281780I 3 * 5 + 248119176s 2 * 6 + 167349456s* 7 + 46204020* 8 ), 
a 7 = -2(158802s 7 + 2143917s 6 * + 15477603s 5 * 2 + 67698320s 4 * 3 + 180289870s 3 * 4 

+280325328s 2 * 5 + 228221280s* 6 + 73941984* 7 ), 
a 8 = 775181s 6 + 9628777s 5 * + 60464369s 4 * 2 + 217547464s 3 * 3 + 438545220s 2 * 4 

+444319344I* 5 + 172576656* 6 , 
a 9 = -2(674202s 5 + 7783209s 4 * + 41993932s 3 * 2 + 117212424s 2 * 3 + 154359000s* 4 

+73984752* 5 ), 
a w = 4(446021s 4 + 4708219s 3 * + 20480415s 2 * 2 + 37508616s* 3 + 23128740* 4 ), 
ail = -8(233734s 3 + 2111409s 2 * + 6071274s* 2 + 5141880* 3 ), 

012 = 6(259913s 2 + 1570956s* + 2057724* 2 ), 

13 = -72(115375 + 31194*), 
oi4 = 187272. 

Note that the long-distance factor (Og ( 3 S\)} 2 does not appear in eq. (A. 2) because it has 
been factored out in eq. (2.1). 

B Parton-level cross sections for gg — y cc ni ( 3 Si) + bb n2 ( 3 Si) 

In this appendix, we present the parton-level differential cross sections for pp — y Hi+H^+X 
for H\ = J/ip and H2 = T. In pp — > J/ip + T + X we consider the three channels 
gg -> cc 8 ( 3 5i) + bb 8 ( 3 S x ), 99 -> cci( 3 Si) + t* 8 ( 3 Si), and gg -* ccs^S^ + 661 ( 3 5i). The 
Mandelstam variables are defined in eq. (A.l), where p\ and p 2 are momenta for cc ni and 
bb n2 , respectively. 
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B.l gg^ccs^S^ + bbsCSi) 

The parton- level differential cross section for gg — > ccs( 3 5'i) + bbs^Si) is given by 

d" 8 

-^[99 -)• ccs^S!) + 6& 8 ( 3 Si)] = F 1 J2 b ijm f /t m^, (B.l) 

i,j=0 



where 



7? a 4 



m? T , , — m 2 ,^ 2 ! 2 



r\ = 7 

108m j/V m T« 2 ( t ~ m j/- ! /-) 2 (* ~ W-x) 2 ^ - m j/^) 2 ('" ~ m x) 2 [s 2 - (" l j/^ - " l T 

(B.2 
Here, b%j = bji and non- vanishing elements of 6j,'s are 

600 = 54s 2 ?(s + t) 2 (s 2 + si + ?) 3 , 

6 01 = -54s?(s + i)(s 2 + s£ + ?)(3s 5 + 6sH + 6s 3 ? + 2s 2 ? - s? - 2?), 

602 = 2? (46s 8 + 49s 7 i - 199s 6 ? - 701s 5 ? - 1090s 4 ? - 1004s 3 ? - 548s 2 ? - 135s? 

+27?), 

603 = 2?(97s 7 + 499s 6 £ + 1147s 5 * 2 + 1520s 4 ? + 1218s 3 ? + 508s 2 ? - 16s? - 135?), 

6 04 = -2?(154s 6 + 484s 5 £ + 636s 4 * 2 + 323s 3 ? - 177s 2 ? - 410s? - 289?), 

6 05 = 2?(49s 5 - 26s 4 t - 350s 3 ? - 634s 2 ? - 613s? - 346?), 

6 6 = 2?(62s 4 + 281s 3 i + 465s 2 ? + 430s? + 249?), 
b 07 = -2? (65s 3 + 149s 2 £ + 1495? + 103?), 

6 8 = 38?(s 2 + st + ?), 

b u = t(54s 9 + 724s 8 t + 2167s 7 ? + 3438s 6 ? + 3305s 5 ? + 1685s 4 ? + 68s 3 ? - 361s 2 ? 

+216* 9 ), 
b 12 = -t(152s 8 + 1025s 7 i + 2436s 6 ? + 3632s 5 ? + 3862s 4 ? + 3159s 3 ? + 2468s 2 ? 

+2074s? + 1350?), 
b 13 = t(34s 7 + 212s 6 t + 1089s 5 ? + 3503s 4 ? + 6172s 3 ? + 7093s 2 ? + 5980s? + 3439?), 
614 = t(200s 6 + 314s 5 £ - 1478s 4 ? - 5353s 3 ? - 8018s 2 ? - 7667s? - 4728?), 
6 1B = -t(18s 5 - 740s 4 i - 3227s 3 ? - 5442s 2 ? - 5707s? - 3962?), 

6 16 = -t(340s 4 + 1469s 3 i + 2662s 2 ? + 2833s? + 2194?), 

617 = t(298s 3 + 780s 2 £+ 893s? + 807?), 

6 18 = -38t(2s 2 + 3s£ + 4?), 

6 22 = 152s 8 + 1243s 7 t + 5142s 6 ? + 12412s 5 ? + 20633s 4 ? + 24264s 3 ? + 20866s 2 ? 

+13616s? + 6546?, 

6 23 = -304s 7 - 2502s 6 £ - 9334s 5 ? - 22006s 4 ? - 34432s 3 ? - 36842s 2 ? - 27395s? 



-14020?, 



624 = 146s 6 + 1808s 5 t + 7883s 4 ? + 19046s 3 ? + 27583s 2 ? + 25828s? + 15958?, 

6 25 = -2(40s 5 + 380s 4 t + 2229s 3 ? + 4889s 2 ? + 6304s? + 5163?), 
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&26 = 216s 4 + 915s 3 t + 2363s 2 £ 2 + 3680s£ 3 + 4078? , 
b 27 = -168s 3 - 666s 2 £ - 891st 2 - 1168? , 
6 28 = 38(s 2 + 3si + 6?), 

6 33 = 1324s 6 + 7637s 5 £ + 23669s 4 ? + 45278s 3 ? + 57940s 2 ? + 48610s? + 27204?, 

634 = -1306s 5 - 8050s 4 t - 22644s 3 ? - 38898s 2 ? - 40961st 4 - 28098?, 

635 = 148s 4 + 3198s 3 i + 9967s 2 ? + 16183s? + 15507?, 

6 36 = -2(4s 3 + 230s 2 £+ 1321st 2 + 2166?), 

637 = 184s 2 + 295si + 722?, 

6 38 = -38(s + 4t), 

644 = 2651s 4 + 10864s 3 £ + 25178s 2 ? + 32284s? + 27180?, 

6 45 = -1267s 3 - 6028s 2 t - 11581st 2 - 13788?, 

6 4 6 = -171s 2 + 1276s* + 2998?, 

647 = s — 138£, 

648 = 38, 

6 55 = 1665s 2 + 3866s* + 6684?, 

6 56 = -341s - 1330? 

657 = -17, 
6 66 = 282. 

B.2 gg^cc 1 ( 3 S 1 ) + bb 8 ( 3 S 1 ) 

The parton-level differential cross section for gg — > cci( 3 Si) + bbs^Sx) is given by 

r 3 

-Jb 5 -* cc^Sx) + bb 8 ( 3 Sx)] = F 2 £ Ci jm %v4, (B.3) 

i,j=0 



where 



107r 3 a 4 



243m J/ ^m 3 f s 2 (t - m 2 j /%jj ) 2 (u - m 2 /v ,) 2 (s - m 2 /v , + m\) 2 



(B.4) 



and non-vanishing elements of 



C{ j s are 



cq! = ?(s + i) 2 , (B.5a) 

c 02 = -2t(s + t) 2 , (B.5b) 

c 03 =t(2s + t), (B.5c) 

cio = 2(s 2 + si + ?) 2 , (B.5d) 

en = _2?(s + 3t), (B.5e) 

C12 = 2(s 2 + 3?), (B.5f) 

c 13 = -2(s + £), (B.5g) 

C20 = -2(s + £)(2s 2 + si + 2?), (B.5h) 
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da 


cc 8 ( 3 S 


i) + bb 1 ( i S 1 )} 


i,j=0 


i m j/tp m y > 


where 


c ji 


's are 


defined in eq. 


(B.5) 


and 
















F 3 








10n 3 a 4 s 








243m 3 7/ ,mrs 2 


(*- 


my) 2 (u - 


- m 2 c ) 2 (s + 


m h - 


- my) 2 



C21 = 3s 2 + 2si + 9P, (B.5i) 

c 22 = 2(2s - 3t), (B.5j) 

C23 = 1, (B.5k) 

c 30 = 2(s 2 + si + P), (B.51) 

c 3 i = -2(s + 2t), (B.5m) 

C32 = 2. (B.5n) 

B.3 M ^cc 8 ( 3 5i) + 66i( 3 5 1 ) 

The parton-level differential cross section for gg — > ccg( 3 Si) + 661 ( 3 «Si) is given by 

(B.6) 



(B.7) 
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